textbook can be found: https://otexts.org/fpp2/
retaildata <- read.csv("https://robjhyndman.com/nyc2018/retail.csv")
mytimeseries <- ts(retaildata[,4], frequency=12, start=c(1982,4))
autoplot(a10)
ggseasonplot(a10)
ggseasonplot(a10, polar=TRUE) + ylab("$ million")
ggsubseriesplot(a10)
beer <- window(ausbeer, start=1992)
autoplot(beer)
ggseasonplot(beer)
ggsubseriesplot(beer)
autoplot(mytimeseries)
ggseasonplot(mytimeseries)
ggseasonplot(mytimeseries, polar=TRUE) + ylab("$ million")
ggsubseriesplot(mytimeseries)
Seasonality is fixed length, cyclic data have variables length
gglagplot(beer, lags=9)
ggAcf(beer)
elec2 <- window(elec, start=1980)
ggAcf(elec2, lag.max=48)
autoplot(goog)
ggAcf(goog, lag.max=100)
wn <- ts(rnorm(36))
autoplot(wn)
ggAcf(wn)
pigs2 <- window(pigs, start=1990)
autoplot(pigs2)
gglagplot(pigs2)
ggAcf(pigs2)
ggseasonplot(pigs2)
ggsubseriesplot(pigs2)
autoplot(bicoal)
gglagplot(bicoal)
ggAcf(bicoal)
autoplot(chicken)
gglagplot(chicken)
ggAcf(chicken)
autoplot(bricksq)
gglagplot(bricksq)
ggAcf(bricksq)
ggseasonplot(bricksq)
ggsubseriesplot(bricksq)
meanf(goog, h=20)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1001 599.4252 446.0804 752.77 364.7754 834.075
## 1002 599.4252 446.0804 752.77 364.7754 834.075
## 1003 599.4252 446.0804 752.77 364.7754 834.075
## 1004 599.4252 446.0804 752.77 364.7754 834.075
## 1005 599.4252 446.0804 752.77 364.7754 834.075
## 1006 599.4252 446.0804 752.77 364.7754 834.075
## 1007 599.4252 446.0804 752.77 364.7754 834.075
## 1008 599.4252 446.0804 752.77 364.7754 834.075
## 1009 599.4252 446.0804 752.77 364.7754 834.075
## 1010 599.4252 446.0804 752.77 364.7754 834.075
## 1011 599.4252 446.0804 752.77 364.7754 834.075
## 1012 599.4252 446.0804 752.77 364.7754 834.075
## 1013 599.4252 446.0804 752.77 364.7754 834.075
## 1014 599.4252 446.0804 752.77 364.7754 834.075
## 1015 599.4252 446.0804 752.77 364.7754 834.075
## 1016 599.4252 446.0804 752.77 364.7754 834.075
## 1017 599.4252 446.0804 752.77 364.7754 834.075
## 1018 599.4252 446.0804 752.77 364.7754 834.075
## 1019 599.4252 446.0804 752.77 364.7754 834.075
## 1020 599.4252 446.0804 752.77 364.7754 834.075
naive(goog, h=20)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1001 813.67 802.4765 824.8634 796.5511 830.7889
## 1002 813.67 797.8401 829.4999 789.4602 837.8798
## 1003 813.67 794.2824 833.0576 784.0192 843.3208
## 1004 813.67 791.2831 836.0569 779.4322 847.9078
## 1005 813.67 788.6407 838.6993 775.3910 851.9490
## 1006 813.67 786.2517 841.0882 771.7374 855.6026
## 1007 813.67 784.0549 843.2851 768.3776 858.9623
## 1008 813.67 782.0101 845.3298 765.2504 862.0896
## 1009 813.67 780.0896 847.2503 762.3133 865.0267
## 1010 813.67 778.2732 849.0668 759.5353 867.8047
## 1011 813.67 776.5455 850.7945 756.8930 870.4470
## 1012 813.67 774.8947 852.4452 754.3684 872.9716
## 1013 813.67 773.3114 854.0285 751.9469 875.3931
## 1014 813.67 771.7879 855.5520 749.6169 877.7231
## 1015 813.67 770.3179 857.0220 747.3688 879.9712
## 1016 813.67 768.8962 858.4438 745.1944 882.1456
## 1017 813.67 767.5182 859.8218 743.0869 884.2530
## 1018 813.67 766.1802 861.1598 741.0406 886.2993
## 1019 813.67 764.8789 862.4611 739.0504 888.2895
## 1020 813.67 763.6114 863.7286 737.1119 890.2280
snaive(goog, h=20)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1001 813.67 802.4765 824.8634 796.5511 830.7889
## 1002 813.67 797.8401 829.4999 789.4602 837.8798
## 1003 813.67 794.2824 833.0576 784.0192 843.3208
## 1004 813.67 791.2831 836.0569 779.4322 847.9078
## 1005 813.67 788.6407 838.6993 775.3910 851.9490
## 1006 813.67 786.2517 841.0882 771.7374 855.6026
## 1007 813.67 784.0549 843.2851 768.3776 858.9623
## 1008 813.67 782.0101 845.3298 765.2504 862.0896
## 1009 813.67 780.0896 847.2503 762.3133 865.0267
## 1010 813.67 778.2732 849.0668 759.5353 867.8047
## 1011 813.67 776.5455 850.7945 756.8930 870.4470
## 1012 813.67 774.8947 852.4452 754.3684 872.9716
## 1013 813.67 773.3114 854.0285 751.9469 875.3931
## 1014 813.67 771.7879 855.5520 749.6169 877.7231
## 1015 813.67 770.3179 857.0220 747.3688 879.9712
## 1016 813.67 768.8962 858.4438 745.1944 882.1456
## 1017 813.67 767.5182 859.8218 743.0869 884.2530
## 1018 813.67 766.1802 861.1598 741.0406 886.2993
## 1019 813.67 764.8789 862.4611 739.0504 888.2895
## 1020 813.67 763.6114 863.7286 737.1119 890.2280
rwf(goog, drift=TRUE, h=50)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1001 814.0912 802.8996 825.2829 796.9751 831.2073
## 1002 814.5125 798.6773 830.3477 790.2946 838.7304
## 1003 814.9338 795.5300 834.3376 785.2582 844.6093
## 1004 815.3550 792.9383 837.7718 781.0716 849.6385
## 1005 815.7763 790.7011 840.8514 777.4271 854.1254
## 1006 816.1976 788.7154 843.6797 774.1673 858.2278
## 1007 816.6188 786.9200 846.3176 771.1984 862.0393
## 1008 817.0401 785.2749 848.8052 768.4595 865.6207
## 1009 817.4613 783.7526 851.1701 765.9083 869.0144
## 1010 817.8826 782.3329 853.4323 763.5140 872.2512
## 1011 818.3039 781.0005 855.6072 761.2533 875.3544
## 1012 818.7251 779.7438 857.7064 759.1083 878.3419
## 1013 819.1464 778.5533 859.7395 757.0646 881.2281
## 1014 819.5676 777.4214 861.7139 755.1106 884.0247
## 1015 819.9889 776.3419 863.6359 753.2366 886.7412
## 1016 820.4102 775.3095 865.5108 751.4347 889.3856
## 1017 820.8314 774.3199 867.3430 749.6982 891.9647
## 1018 821.2527 773.3692 869.1362 748.0212 894.4842
## 1019 821.6739 772.4542 870.8937 746.3988 896.9491
## 1020 822.0952 771.5720 872.6184 744.8267 899.3638
## 1021 822.5165 770.7202 874.3127 743.3010 901.7320
## 1022 822.9377 769.8966 875.9788 741.8184 904.0571
## 1023 823.3590 769.0993 877.6187 740.3759 906.3420
## 1024 823.7803 768.3265 879.2340 738.9710 908.5895
## 1025 824.2015 767.5766 880.8264 737.6012 910.8019
## 1026 824.6228 766.8483 882.3973 736.2643 912.9812
## 1027 825.0440 766.1403 883.9478 734.9586 915.1295
## 1028 825.4653 765.4515 885.4791 733.6821 917.2485
## 1029 825.8866 764.7808 886.9923 732.4333 919.3398
## 1030 826.3078 764.1272 888.4884 731.2108 921.4048
## 1031 826.7291 763.4900 889.9682 730.0132 923.4450
## 1032 827.1503 762.8682 891.4325 728.8393 925.4614
## 1033 827.5716 762.2611 892.8821 727.6879 927.4553
## 1034 827.9929 761.6682 894.3176 726.5580 929.4278
## 1035 828.4141 761.0886 895.7397 725.4486 931.3797
## 1036 828.8354 760.5218 897.1489 724.3588 933.3119
## 1037 829.2566 759.9674 898.5459 723.2879 935.2254
## 1038 829.6779 759.4247 899.9311 722.2349 937.1209
## 1039 830.0992 758.8933 901.3050 721.1992 938.9991
## 1040 830.5204 758.3728 902.6681 720.1801 940.8608
## 1041 830.9417 757.8626 904.0208 719.1769 942.7065
## 1042 831.3630 757.3625 905.3634 718.1891 944.5368
## 1043 831.7842 756.8721 906.6963 717.2160 946.3524
## 1044 832.2055 756.3910 908.0200 716.2572 948.1537
## 1045 832.6267 755.9188 909.3346 715.3121 949.9413
## 1046 833.0480 755.4554 910.6406 714.3803 951.7157
## 1047 833.4693 755.0003 911.9382 713.4613 953.4772
## 1048 833.8905 754.5533 913.2277 712.5547 955.2263
## 1049 834.3118 754.1142 914.5094 711.6601 956.9634
## 1050 834.7330 753.6826 915.7835 710.7771 958.6890
goog %>% meanf(h=20) %>% autoplot
goog %>% naive(h=20) %>% autoplot
goog %>% snaive(h=20) %>% autoplot
goog %>% rwf(drift=TRUE, h=50) %>% autoplot
If the fitted values aren’t good, then
goog %>% rwf(drift=TRUE, h=500) %>% fitted() -> z
autoplot(goog, series='Data') + autolayer(z, series='Fitted')
## Warning: Removed 1 rows containing missing values (geom_path).
Residuals \(e_t\) should not be correlated, if they are then information left in the residuals shoudl be used in the forecast
If the mean is not 0 then the forecast is biased
fits <- fitted(naive(goog200))
autoplot(goog200, series='Data') + autolayer(fits, series='Fitted')
## Warning: Removed 1 rows containing missing values (geom_path).
res <- residuals(naive(goog200))
autoplot(res)
gghistogram(res, add.normal=TRUE)
## Warning: Removed 1 rows containing non-finite values (stat_bin).
ggAcf(res)
checkresiduals(naive(goog200))
##
## Ljung-Box test
##
## data: Residuals from Naive method
## Q* = 11.031, df = 10, p-value = 0.3551
##
## Model df: 0. Total lags used: 10
checkresiduals also computes the Ljung-Box test to make sure the residuals are white noise, if p-value is large it is white noise
beer <- window(ausbeer, start=1992)
fc <- snaive(beer)
autoplot(fc)
res <- residuals(fc)
autoplot(res)
checkresiduals(fc)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 32.269, df = 8, p-value = 8.336e-05
##
## Model df: 0. Total lags used: 8
This tiny p-value means it’s not white noise and we can do better.
\[ MAPE = 100mean(|e_{T+h}|/y_{T+h}) \]
MAPE doesn’t work well when numbers are close to 0 because of division by \(y_{T+h}\)
MAPE also doesn’t make sense when there is no natura 0, like temperature
Instead use Mean Absolute Scaled Error
\[ MASE = T^{-1} \sum_{t=1}^T |y_t \hat{y_{t|t-1}}| / Q \]
where Q is a stable measure of the scale of the time series
googtrain <- window(goog200, end=180)
googfc1 <- meanf(googtrain, h=20)
googfc2 <- naive(googtrain, h=20)
googfc3 <- rwf(googtrain, drift=TRUE, h=20)
accuracy(googfc1, goog200)
## ME RMSE MAE MPE MAPE
## Training set -1.360278e-14 28.74910 20.06324 -0.4195606 4.582412
## Test set 8.242659e+01 82.89387 82.42659 15.9262898 15.926290
## MASE ACF1 Theil's U
## Training set 5.261155 0.9547241 NA
## Test set 21.614609 0.7822826 20.39543
accuracy(googfc2, goog200)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.602728 6.404477 3.813467 0.1256394 0.8712766 1.000000
## Test set 16.041900 18.291917 16.041900 3.0762438 3.0762438 4.206645
## ACF1 Theil's U
## Training set -0.06213755 NA
## Test set 0.78228260 4.539165
accuracy(googfc3, goog200)
## ME RMSE MAE MPE MAPE
## Training set -2.191158e-14 6.376052 3.887099 -0.01363486 0.8882792
## Test set 9.713257e+00 11.335822 9.713257 1.86155541 1.8615554
## MASE ACF1 Theil's U
## Training set 1.019309 -0.06213755 NA
## Test set 2.547094 0.70229763 2.811493
train <- window(mytimeseries, end=c(2010, 12))
test <- window(mytimeseries, start=2011)
autoplot(cbind(Training=train, Test=test))
fcst1 <- snaive(train, h=length(test))
accuracy(fcst1, test)
## ME RMSE MAE MPE MAPE MASE
## Training set 5.502402 27.74935 19.33784 3.369450 10.447161 1.0000000
## Test set 5.227586 23.40458 18.60000 1.619239 8.172045 0.9618449
## ACF1 Theil's U
## Training set 0.8703252 NA
## Test set 0.4544630 0.9518744
checkresiduals(fcst1)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 1256.8, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24
fcst1 %>% autoplot() + autolayer(test)
Rob coined the phrase
rolling training data, all starting at the same point
e <- tsCV(goog200, rwf, drift=TRUE, h=1)
sqrt(mean(e^2, na.rm=TRUE))
## [1] 6.233245
Now with pipes
e <- goog200 %>% tsCV(forecastfunction=rwf, drift=TRUE, h=1)
e^2 %>% mean(na.rm=TRUE) %>% sqrt()
## [1] 6.233245
Now get error for h=1, h=2, h=3, etc
e <- goog200 %>% tsCV(forecastfunction=rwf, drift=TRUE, h=12)
e^2 %>% colMeans(na.rm=TRUE) %>% sqrt()
## h=1 h=2 h=3 h=4 h=5 h=6 h=7
## 6.233245 8.634156 10.833826 13.032527 14.962450 16.483554 18.108701
## h=8 h=9 h=10 h=11 h=12
## 19.949726 21.648751 23.237067 24.635515 25.870110
e1 <- mytimeseries %>% tsCV(forecastfunction=rwf, drift=TRUE, h=12)
e2 <- mytimeseries %>% tsCV(forecastfunction=meanf, h=12)
e3 <- mytimeseries %>% tsCV(forecastfunction=naive, h=12)
e4 <- mytimeseries %>% tsCV(forecastfunction=snaive, h=12)
e1^2 %>% colMeans(na.rm=TRUE) %>% sqrt()
## h=1 h=2 h=3 h=4 h=5 h=6 h=7 h=8
## 20.04797 23.06671 24.18534 25.51110 27.27162 30.09422 29.46610 29.93850
## h=9 h=10 h=11 h=12
## 30.69911 31.68924 31.27128 26.49646
e2^2 %>% colMeans(na.rm=TRUE) %>% sqrt()
## h=1 h=2 h=3 h=4 h=5 h=6 h=7 h=8
## 67.49756 67.84488 68.19028 68.53594 68.88044 69.22489 69.56606 69.91056
## h=9 h=10 h=11 h=12
## 70.25265 70.59840 70.94401 71.29212
e3^2 %>% colMeans(na.rm=TRUE) %>% sqrt()
## h=1 h=2 h=3 h=4 h=5 h=6 h=7 h=8
## 20.01105 22.98583 24.06746 25.36152 27.08889 29.83772 29.20396 29.64566
## h=9 h=10 h=11 h=12
## 30.36709 31.31045 30.89591 26.30046
e4^2 %>% colMeans(na.rm=TRUE) %>% sqrt()
## h=1 h=2 h=3 h=4 h=5 h=6 h=7 h=8
## 26.30046 26.33190 26.36324 26.39421 26.42471 26.45617 26.48818 26.52014
## h=9 h=10 h=11 h=12
## 26.55237 26.58396 26.61484 26.64629
For one-ahead, this gives the overall error for each method
e1 <- mytimeseries %>% tsCV(forecastfunction=rwf, drift=TRUE, h=1)
e2 <- mytimeseries %>% tsCV(forecastfunction=meanf, h=1)
e3 <- mytimeseries %>% tsCV(forecastfunction=naive, h=1)
e4 <- mytimeseries %>% tsCV(forecastfunction=snaive, h=1)
e1^2 %>% mean(na.rm=TRUE) %>% sqrt()
## [1] 20.04797
e2^2 %>% mean(na.rm=TRUE) %>% sqrt()
## [1] 67.49756
e3^2 %>% mean(na.rm=TRUE) %>% sqrt()
## [1] 20.01105
e4^2 %>% mean(na.rm=TRUE) %>% sqrt()
## [1] 26.30046
cbind(RWF=e1^2, Mean=e2^2, Naive=e3^2, SNaive=e4^2) %>% na.omit() %>% colMeans()
## RWF Mean Naive SNaive
## 411.4290 4685.8665 410.0438 691.7144
Robert Goodall Brown
fc <- ses(oil, h=12)
autoplot(fc)
summary(fc)
##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = oil, h = 12)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 110.8832
##
## sigma: 49.0507
##
## AIC AICc BIC
## 576.1569 576.6903 581.8324
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 8.806147 48.03921 35.52053 2.135668 10.73828 0.9796835
## ACF1
## Training set 0.1646214
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2014 542.3412 479.4802 605.2022 446.2037 638.4788
## 2015 542.3412 453.4468 631.2356 406.3890 678.2935
## 2016 542.3412 433.4701 651.2124 375.8372 708.8453
## 2017 542.3412 416.6287 668.0537 350.0805 734.6019
## 2018 542.3412 401.7911 682.8914 327.3883 757.2941
## 2019 542.3412 388.3767 696.3057 306.8729 777.8096
## 2020 542.3412 376.0410 708.6415 288.0069 796.6755
## 2021 542.3412 364.5591 720.1233 270.4469 814.2355
## 2022 542.3412 353.7751 730.9074 253.9542 830.7283
## 2023 542.3412 343.5753 741.1072 238.3549 846.3275
## 2024 542.3412 333.8739 750.8085 223.5180 861.1644
## 2025 542.3412 324.6044 760.0781 209.3415 875.3410
summary(fc$model)
## Simple exponential smoothing
##
## Call:
## ses(y = oil, h = 12)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 110.8832
##
## sigma: 49.0507
##
## AIC AICc BIC
## 576.1569 576.6903 581.8324
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 8.806147 48.03921 35.52053 2.135668 10.73828 0.9796835
## ACF1
## Training set 0.1646214
Charles Holt (US Navy) (1957)
window(ausair, start=1990, end=2004) %>%
holt(h=5, PI=FALSE) %>%
summary()
##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = ., h = 5, PI = FALSE)
##
## Smoothing parameters:
## alpha = 1e-04
## beta = 1e-04
##
## Initial states:
## l = 19.1015
## b = 1.4307
##
## sigma: 1.9943
##
## AIC AICc BIC
## 66.67669 73.34336 70.21694
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.09267432 1.707786 1.474578 -0.8251858 5.100699 0.8045133
## ACF1
## Training set 0.4475061
##
## Forecasts:
## Point Forecast
## 2005 41.99044
## 2006 43.42100
## 2007 44.85155
## 2008 46.28211
## 2009 47.71267
window(ausair, start=1990, end=2004) %>%
holt(h=5, PI=FALSE) %>%
autoplot()
1980s: Gardner (Houston) & McKinsey (Scotland)
autoplot(livestock)
livestock2 <- livestock %>% window(start=1970, end=2000)
livestock2 %>% holt(h=20, damped=TRUE) %>% autoplot()
Beta=0.3 is pretty big
phi=0.8 is artifically constrained to be the smallest phi can be
SES for non-trended, Holt for trended
eggs2 <- window(eggs, start=1900, end=1993)
fc1 <- ses(eggs2, h=100)
fc2 <- holt(eggs2, h=100)
fc3 <- holt(eggs2, h=100, damped=TRUE)
fc1 %>% autoplot()
fc2 %>% autoplot()
fc3 %>% autoplot()
fc1 %>% autoplot(PI=FALSE)
fc2 %>% autoplot(PI=FALSE)
fc3 %>% autoplot(PI=FALSE)
residuals(fc1)^2 %>% mean() %>% sqrt()
## [1] 26.56395
residuals(fc2)^2 %>% mean() %>% sqrt()
## [1] 26.58219
residuals(fc3)^2 %>% mean() %>% sqrt()
## [1] 26.54019
fc1 %>% accuracy()
## ME RMSE MAE MPE MAPE MASE
## Training set -2.740622 26.56395 19.38206 -2.872958 10.05648 0.9560879
## ACF1
## Training set -0.005983602
fc2 %>% accuracy()
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
## ACF1
## Training set 0.01348202
fc3 %>% accuracy()
## ME RMSE MAE MPE MAPE MASE
## Training set -2.891496 26.54019 19.2795 -2.907633 10.01894 0.9510287
## ACF1
## Training set -0.003195358
list(SES=fc1, Trend=fc2, DampedTrend=fc3) %>% purrr::map_df(~ accuracy(.x) %>% as.data.frame(), .id='Model')
## Model ME RMSE MAE MPE MAPE MASE
## 1 SES -2.74062165 26.56395 19.38206 -2.872958 10.056482 0.9560879
## 2 Trend 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
## 3 DampedTrend -2.89149552 26.54019 19.27950 -2.907633 10.018944 0.9510287
## ACF1
## 1 -0.005983602
## 2 0.013482023
## 3 -0.003195358
fc1 %>% checkresiduals()
##
## Ljung-Box test
##
## data: Residuals from Simple exponential smoothing
## Q* = 10.475, df = 8, p-value = 0.2333
##
## Model df: 2. Total lags used: 10
fc2 %>% checkresiduals()
##
## Ljung-Box test
##
## data: Residuals from Holt's method
## Q* = 11.026, df = 6, p-value = 0.08759
##
## Model df: 4. Total lags used: 10
fc3 %>% checkresiduals()
##
## Ljung-Box test
##
## data: Residuals from Damped Holt's method
## Q* = 10.782, df = 5, p-value = 0.05587
##
## Model df: 5. Total lags used: 10
Holt Winters: Holt wrote the pape in 1957, separately Winters (Holt’s students) wrote the code 3 years later
aust <- window(austourists, start=2005)
fc1 <- hw(aust, seasonal='additive')
fc2 <- hw(aust, seasonal='multiplicative')
fc1 %>% autoplot()
fc2 %>% autoplot()
fc1 %>% summary()
##
## Forecast method: Holt-Winters' additive method
##
## Model Information:
## Holt-Winters' additive method
##
## Call:
## hw(y = aust, seasonal = "additive")
##
## Smoothing parameters:
## alpha = 0.3063
## beta = 1e-04
## gamma = 0.4263
##
## Initial states:
## l = 32.2597
## b = 0.7014
## s=1.3106 -1.6935 -9.3132 9.6962
##
## sigma: 1.9494
##
## AIC AICc BIC
## 234.4171 239.7112 250.4748
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.008115785 1.763305 1.374062 -0.2860248 2.973922 0.4502579
## ACF1
## Training set -0.06272507
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2016 Q1 76.09837 73.60011 78.59664 72.27761 79.91914
## 2016 Q2 51.60333 48.99039 54.21626 47.60718 55.59947
## 2016 Q3 63.96867 61.24582 66.69153 59.80443 68.13292
## 2016 Q4 68.37170 65.54313 71.20027 64.04578 72.69762
## 2017 Q1 78.90404 75.53440 82.27369 73.75061 84.05747
## 2017 Q2 54.40899 50.95325 57.86473 49.12389 59.69409
## 2017 Q3 66.77434 63.23454 70.31414 61.36069 72.18799
## 2017 Q4 71.17737 67.55541 74.79933 65.63806 76.71667
fc2 %>% summary()
##
## Forecast method: Holt-Winters' multiplicative method
##
## Model Information:
## Holt-Winters' multiplicative method
##
## Call:
## hw(y = aust, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha = 0.4406
## beta = 0.0134
## gamma = 0.0023
##
## Initial states:
## l = 32.4875
## b = 0.6974
## s=1.0237 0.9618 0.7704 1.2442
##
## sigma: 0.0367
##
## AIC AICc BIC
## 221.1313 226.4254 237.1890
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.09206228 1.575631 1.25496 -0.0006505533 2.70539 0.4112302
## ACF1
## Training set -0.07955726
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2016 Q1 80.08894 76.31865 83.85922 74.32278 85.85509
## 2016 Q2 50.15482 47.56655 52.74309 46.19640 54.11324
## 2016 Q3 63.34322 59.80143 66.88502 57.92652 68.75993
## 2016 Q4 68.17810 64.08399 72.27221 61.91670 74.43950
## 2017 Q1 83.80112 78.43079 89.17146 75.58790 92.01434
## 2017 Q2 52.45291 48.88795 56.01787 47.00077 57.90504
## 2017 Q3 66.21274 61.46194 70.96353 58.94702 73.47845
## 2017 Q4 71.23205 65.85721 76.60690 63.01194 79.45217
fc1 %>% checkresiduals()
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' additive method
## Q* = 8.6117, df = 3, p-value = 0.03493
##
## Model df: 8. Total lags used: 11
fc2 %>% checkresiduals()
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 6.1188, df = 3, p-value = 0.106
##
## Model df: 8. Total lags used: 11
fc1$model %>% autoplot()
fc2$model %>% autoplot()
F Gardner uses mostly this
So does Rob if he can only use 1 tool
If you only have 1 tool this is it
fc1 <- mytimeseries %>% hw(h=12, seasonal='multiplicative', damped=FALSE)
fc2 <- mytimeseries %>% hw(h=12, seasonal='multiplicative', damped=TRUE)
fc1 %>% autoplot()
fc2 %>% autoplot()
fc1$model %>% autoplot()
fc2$model %>% autoplot()
fc1 %>% summary()
##
## Forecast method: Holt-Winters' multiplicative method
##
## Model Information:
## Holt-Winters' multiplicative method
##
## Call:
## hw(y = ., h = 12, seasonal = "multiplicative", damped = FALSE)
##
## Smoothing parameters:
## alpha = 0.6979
## beta = 0.0087
## gamma = 1e-04
##
## Initial states:
## l = 63.7459
## b = 1.0282
## s=0.9921 0.933 0.9936 1.2091 1.013 1.0172
## 0.979 0.995 0.9791 0.9381 0.9746 0.9761
##
## sigma: 0.0468
##
## AIC AICc BIC
## 4414.235 4415.713 4483.398
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.2054187 9.280748 6.244629 -0.2790438 3.338901 0.3336358
## ACF1
## Training set 0.05182697
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Apr 2018 245.7755 231.0502 260.5008 223.2551 268.2959
## May 2018 245.7272 227.7016 263.7529 218.1593 273.2952
## Jun 2018 236.8256 216.7102 256.9411 206.0617 267.5896
## Jul 2018 247.4718 223.8776 271.0659 211.3876 283.5559
## Aug 2018 251.7887 225.3713 278.2060 211.3868 292.1905
## Sep 2018 248.0573 219.8081 276.3064 204.8539 291.2606
## Oct 2018 258.0681 226.4893 289.6470 209.7724 306.3638
## Nov 2018 257.3264 223.7535 290.8993 205.9811 308.6717
## Dec 2018 307.4929 264.9803 350.0056 242.4754 372.5105
## Jan 2019 253.0461 216.1567 289.9354 196.6287 309.4634
## Feb 2019 237.8969 201.4787 274.3151 182.2000 293.5938
## Mar 2019 253.2718 212.6985 293.8451 191.2203 315.3233
fc2 %>% summary()
##
## Forecast method: Damped Holt-Winters' multiplicative method
##
## Model Information:
## Damped Holt-Winters' multiplicative method
##
## Call:
## hw(y = ., h = 12, seasonal = "multiplicative", damped = TRUE)
##
## Smoothing parameters:
## alpha = 0.7888
## beta = 1e-04
## gamma = 1e-04
## phi = 0.9738
##
## Initial states:
## l = 63.9173
## b = 0.2053
## s=0.9948 0.935 0.9961 1.2126 1.0159 1.0202
## 0.9781 0.992 0.9751 0.9339 0.9713 0.975
##
## sigma: 0.0472
##
## AIC AICc BIC
## 4418.242 4419.898 4491.474
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.5050393 9.194427 6.128123 0.2458937 3.290157 0.3274111
## ACF1
## Training set -0.04640084
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Apr 2018 245.8588 230.9987 260.7189 223.1323 268.5854
## May 2018 244.9277 226.0636 263.7917 216.0775 273.7778
## Jun 2018 235.5175 214.1719 256.8632 202.8722 268.1629
## Jul 2018 245.8983 220.7007 271.0959 207.3620 284.4347
## Aug 2018 250.1604 221.8689 278.4519 206.8922 293.4286
## Sep 2018 246.6482 216.3561 276.9404 200.3204 292.9760
## Oct 2018 257.2616 223.3454 291.1778 205.3912 309.1320
## Nov 2018 256.1870 220.2465 292.1275 201.2208 311.1532
## Dec 2018 305.7881 260.4472 351.1290 236.4452 375.1310
## Jan 2019 251.2252 212.0682 290.3821 191.3398 311.1106
## Feb 2019 235.7881 197.3285 274.2477 176.9692 294.6070
## Mar 2019 250.8599 208.1987 293.5212 185.6152 316.1047
fc1 %>% checkresiduals()
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 41.483, df = 8, p-value = 1.693e-06
##
## Model df: 16. Total lags used: 24
fc2 %>% checkresiduals()
##
## Ljung-Box test
##
## data: Residuals from Damped Holt-Winters' multiplicative method
## Q* = 42.336, df = 7, p-value = 4.479e-07
##
## Model df: 17. Total lags used: 24
fc1 %>% accuracy()
## ME RMSE MAE MPE MAPE MASE
## Training set -0.2054187 9.280748 6.244629 -0.2790438 3.338901 0.3336358
## ACF1
## Training set 0.05182697
fc2 %>% accuracy()
## ME RMSE MAE MPE MAPE MASE
## Training set 0.5050393 9.194427 6.128123 0.2458937 3.290157 0.3274111
## ACF1
## Training set -0.04640084
Exponential smoothing methods: * return point forecasts
Innovations state space models: * general pobabalistic stuff
estimating ETS models:
ets(mytimeseries)
## ETS(M,Ad,M)
##
## Call:
## ets(y = mytimeseries)
##
## Smoothing parameters:
## alpha = 0.7718
## beta = 0.0025
## gamma = 0.0561
## phi = 0.98
##
## Initial states:
## l = 64.0899
## b = 0.2813
## s=0.9869 0.9362 1.0043 1.1705 1.0216 1.0293
## 0.9848 0.9998 0.993 0.9436 0.9694 0.9606
##
## sigma: 0.0464
##
## AIC AICc BIC
## 4404.041 4405.697 4477.272
ets(mytimeseries) %>% autoplot()
ets(mytimeseries) %>% forecast() %>% autoplot()
USAccDeaths %>%
ets %>%
forecast ->
forecasted_stuff
head(forecasted_stuff)
## $model
## ETS(A,N,A)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.5946
## gamma = 0.002
##
## Initial states:
## l = 9248.3628
## s=-51.3449 -255.3528 218.2901 -121.771 970.7387 1683.237
## 756.092 306.4212 -489.5627 -739.9004 -1537.792 -739.0552
##
## sigma: 292.6907
##
## AIC AICc BIC
## 1140.145 1148.716 1174.295
##
## $mean
## Jan Feb Mar Apr May Jun Jul
## 1979 8397.497 7599.221 8396.595 8646.510 9443.412 9893.065 10819.124
## 1980 8397.497 7599.221 8396.595 8646.510 9443.412 9893.065 10819.124
## Aug Sep Oct Nov Dec
## 1979 10107.105 9015.233 9354.407 8880.418 9085.492
## 1980 10107.105 9015.233 9354.407 8880.418 9085.492
##
## $level
## [1] 80 95
##
## $x
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov
## 1973 9007 8106 8928 9137 10017 10826 11317 10744 9713 9938 9161
## 1974 7750 6981 8038 8422 8714 9512 10120 9823 8743 9129 8710
## 1975 8162 7306 8124 7870 9387 9556 10093 9620 8285 8466 8160
## 1976 7717 7461 7767 7925 8623 8945 10078 9179 8037 8488 7874
## 1977 7792 6957 7726 8106 8890 9299 10625 9302 8314 8850 8265
## 1978 7836 6892 7791 8192 9115 9434 10484 9827 9110 9070 8633
## Dec
## 1973 8927
## 1974 8680
## 1975 8034
## 1976 8647
## 1977 8796
## 1978 9240
##
## $upper
## 80% 95%
## Jan 1979 8772.595 8971.160
## Feb 1979 8035.616 8266.630
## Mar 1979 8886.679 9146.115
## Apr 1979 9184.957 9469.994
## May 1979 10026.222 10334.743
## Jun 1979 10517.092 10847.432
## Jul 1979 11481.809 11832.614
## Aug 1979 10806.314 11176.454
## Sep 1979 9749.151 10137.664
## Oct 1979 10121.466 10527.522
## Nov 1979 9679.242 10102.115
## Dec 1979 9915.072 10354.225
## Jan 1980 9256.534 9711.281
## Feb 1980 8486.738 8956.562
## Mar 1980 9311.707 9796.138
## Apr 1980 9588.408 10087.018
## May 1980 10411.355 10923.753
## Jun 1980 10886.371 11412.195
## Jul 1980 11837.160 12376.076
## Aug 1980 11149.285 11700.983
## Sep 1980 10081.011 10645.200
## Oct 1980 10443.272 11019.681
## Nov 1980 9991.889 10580.266
## Dec 1980 10219.269 10819.454
##
## $lower
## 80% 95%
## Jan 1979 8022.399 7823.834
## Feb 1979 7162.825 6931.812
## Mar 1979 7906.510 7647.075
## Apr 1979 8108.063 7823.026
## May 1979 8860.602 8552.081
## Jun 1979 9269.038 8938.698
## Jul 1979 10156.438 9805.634
## Aug 1979 9407.895 9037.756
## Sep 1979 8281.314 7892.801
## Oct 1979 8587.349 8181.293
## Nov 1979 8081.593 7658.721
## Dec 1979 8255.912 7816.759
## Jan 1980 7538.459 7083.713
## Feb 1980 6711.703 6241.880
## Mar 1980 7481.483 6997.052
## Apr 1980 7704.612 7206.001
## May 1980 8475.469 7963.070
## Jun 1980 8899.759 8373.935
## Jul 1980 9801.087 9262.171
## Aug 1980 9064.924 8513.227
## Sep 1980 7949.455 7385.266
## Oct 1980 8265.543 7689.133
## Nov 1980 7768.947 7180.570
## Dec 1980 7951.715 7351.530
forecasted_stuff %>%
autoplot
h02 %>% ets %>% forecast %>% autoplot
h02 %>% ets %>% print
## ETS(M,Ad,M)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.1953
## beta = 1e-04
## gamma = 1e-04
## phi = 0.9798
##
## Initial states:
## l = 0.3945
## b = 0.0085
## s=0.874 0.8197 0.7644 0.7693 0.6941 1.2838
## 1.326 1.1765 1.1621 1.0955 1.0422 0.9924
##
## sigma: 0.0676
##
## AIC AICc BIC
## -122.90601 -119.20871 -63.17985
autoplot(a10)
a10 %>% ets %>% autoplot
jdl_model <- function(dataset) {
dataset %>% ets -> fit_ets
print(fit_ets)
print(fit_ets %>% autoplot)
print(fit_ets %>% forecast %>% autoplot)
}
jdl_model(bicoal)
## ETS(M,N,N)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.8205
##
## Initial states:
## l = 542.665
##
## sigma: 0.1262
##
## AIC AICc BIC
## 595.2499 595.7832 600.9253
jdl_model(chicken)
## ETS(M,N,N)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.98
##
## Initial states:
## l = 159.8322
##
## sigma: 0.1691
##
## AIC AICc BIC
## 635.2382 635.6018 641.9836
jdl_model(dole)
## ETS(M,Ad,M)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.697
## beta = 0.124
## gamma = 0.303
## phi = 0.902
##
## Initial states:
## l = 2708.6621
## b = 836.017
## s=1.0404 0.8893 0.9103 1.0301 1.0576 1.0584
## 0.9801 0.9632 1.021 0.9838 1.0145 1.0514
##
## sigma: 0.0935
##
## AIC AICc BIC
## 10602.67 10604.30 10676.19
jdl_model(usdeaths)
## ETS(A,N,A)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.5972
## gamma = 0.0019
##
## Initial states:
## l = 9195.6403
## s=-62.6129 -270.0351 263.3823 -89.4907 1005.529 1662.647
## 795.2585 333.326 -551.161 -737.5102 -1552.872 -796.4611
##
## sigma: 294.4663
##
## AIC AICc BIC
## 1141.016 1149.587 1175.166
jdl_model(bricksq)
## ETS(M,A,M)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.973
## beta = 0.0032
## gamma = 1e-04
##
## Initial states:
## l = 196.0437
## b = 4.9839
## s=1.0041 1.0598 1.0269 0.9093
##
## sigma: 0.0514
##
## AIC AICc BIC
## 1725.473 1726.714 1752.864
jdl_model(lynx)
## ETS(M,N,N)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 2372.8047
##
## sigma: 0.9594
##
## AIC AICc BIC
## 2058.138 2058.356 2066.346
jdl_model(ibmclose)
## ETS(A,N,N)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 459.9339
##
## sigma: 7.2637
##
## AIC AICc BIC
## 3648.450 3648.515 3660.182
jdl_model(eggs)
## ETS(M,N,N)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.8198
##
## Initial states:
## l = 278.8889
##
## sigma: 0.1355
##
## AIC AICc BIC
## 1043.286 1043.553 1050.916
jdl_model(ausbeer)
## ETS(M,A,M)
##
## Call:
## ets(y = .)
##
## Smoothing parameters:
## alpha = 0.2087
## beta = 0.0302
## gamma = 0.1984
##
## Initial states:
## l = 255.892
## b = 0.6795
## s=1.183 0.9113 0.8614 1.0444
##
## sigma: 0.0361
##
## AIC AICc BIC
## 2354.249 2355.114 2384.709
## works but needs to print name on every run of the function
data_list <- list(bicoal, chicken, dole, usdeaths, bricksq, lynx, ibmclose, eggs, ausbeer)
data_list$names <- c("bicoal", "chicken", "dole", "usdeaths", "bricksq", "lynx", "ibmclose", "eggs", "ausbeer")
lapply(data_list, jdl_model)
autoplot(mytimeseries)
train <- window(mytimeseries, end=c(2010,12))
test <- window(mytimeseries, start=2011)
f1 <- snaive(train, h=length(test))
f2 <- rwf(train, h=length(test))
f3 <- rwf(train, drift = TRUE, h=length(test))
f4 <- meanf(train, h=length(test))
f5 <- hw(train, h=length(test), seasonal = 'multiplicative')
f6 <- ets(train) %>% forecast(h=length(test))
c(
SN=accuracy(f1, test)[2,"RMSE"],
rwf1=accuracy(f2, test)[2,"RMSE"],
rwf2=accuracy(f3, test)[2,"RMSE"],
meanf=accuracy(f4, test)[2,"RMSE"],
hw=accuracy(f5, test)[2,"RMSE"],
ets=accuracy(f6, test)[2,"RMSE"]
)
## SN rwf1 rwf2 meanf hw ets
## 23.40458 37.09060 56.64803 58.73867 59.85277 25.43061
now let’s do time series cross validation
e1 <- tsCV(mytimeseries, snaive, h=12)
e2 <- tsCV(mytimeseries, naive, h=12)
e3 <- tsCV(mytimeseries, rwf, drift=TRUE, h=12)
e4 <- tsCV(mytimeseries, meanf, h=12)
e5 <- tsCV(mytimeseries, hw, h=12, seasonal='multiplicative')
etsfc <- function(y,h){
y %>% ets(model="MAM", damped=TRUE) %>%
forecast(h=h)
}
e6 <- tsCV(mytimeseries, etsfc, h=12)
MSE <- cbind(
h=1:12,
SNaive=colMeans(tail(e1, -14)^2, na.rm=TRUE),
Naive=colMeans(tail(e2, -14)^2, na.rm=TRUE),
drift=colMeans(tail(e3, -14)^2, na.rm=TRUE),
Mean=colMeans(tail(e4, -14)^2, na.rm=TRUE),
HW=colMeans(tail(e5, -14)^2, na.rm=TRUE),
ETS=colMeans(tail(e6, -14)^2, na.rm=TRUE)
)
MSE
## h SNaive Naive drift Mean HW ETS
## h=1 1 695.0202 411.9992 413.3894 4708.335 106.7760 104.3003
## h=2 2 696.6544 543.7812 547.3381 4757.235 171.8868 158.4553
## h=3 3 698.2653 596.3238 601.7417 4806.157 251.5070 225.1077
## h=4 4 699.9288 661.9810 668.9494 4855.374 323.9946 294.1775
## h=5 5 701.6235 755.7454 765.1886 4904.751 398.8013 364.3042
## h=6 6 703.3177 917.1529 931.6504 4954.304 479.8465 418.0667
## h=7 7 705.0283 878.4346 893.1175 5003.275 559.1300 494.1943
## h=8 8 706.7069 905.6455 922.2475 5053.316 653.4687 555.8380
## h=9 9 708.3497 950.4590 969.6496 5103.750 730.6287 624.1344
## h=10 10 710.0246 1010.2253 1032.6967 5154.446 822.1079 694.4674
## h=11 11 711.6948 983.8066 1005.4918 5205.465 916.4102 757.7232
## h=12 12 713.3165 713.3165 722.5366 5257.032 996.6152 808.3827
Transformations. Box Cox
elec %>% autoplot
(lambda <- BoxCox.lambda(elec))
## [1] 0.2654076
elec %>%
BoxCox(lambda) %>%
autoplot
fc <- rwf(eggs, drift=TRUE, lambda=0, h=50, level=80)
fc2 <- rwf(eggs, drift=TRUE, lambda=0, h=50, level=80,
biasadj=TRUE)
autoplot(eggs) +
autolayer(fc, series="Simple back transformation") +
autolayer(fc2, series="Bias adjusted", PI=FALSE) +
guides(colour=guide_legend(title="Forecast"))
lambda_adjusted <- function(data) {
data %>% autoplot
lambda <- BoxCox.lambda(data)
print(lambda)
print("with transform")
data %>%
BoxCox(lambda) %>%
forecast %>%
autoplot
print("without transform")
data %>%
forecast %>%
autoplot
}
data_list <- list(usnetelec,
mcopper,
enplanements,
a10
)
lapply(data_list, lambda_adjusted)
## [1] 0.5167714
## [1] "with transform"
## [1] "without transform"
## [1] 0.1919047
## [1] "with transform"
## [1] "without transform"
## [1] -0.2269461
## [1] "with transform"
## [1] "without transform"
## [1] 0.1313326
## [1] "with transform"
## [1] "without transform"
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]